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The diagram technique for the one-band Hubbard model is formulated for the case of moderate to 
strong Hubbard repulsion. The expansion in powers of the hopping constant is expressed in terms 
of site cumulants of electron creation and annihilation operators. For Green's function an equation 
of the Larkin type is derived and solved in a one-loop approximation for the case of two dimensions, 
nearest-neighbor hopping and half-filling. The obtained four-band structure of the spectrum and the 
shapes of the spectral function are close to those observed in Monte Carlo calculations. It is shown 
that the maxima forming the bands are of a dissimilar origin in different regions of the Brillouin 
zone. 

PACS numbers; 71.10.Fd, 71.10.-w 



I 

o 
o 



m 
m 

o 

o 



I 

o 
o 



X 



I. INTRODUCTION 

Systems with strong electron correlations which, in 
particular, cuprate perovskites belong to are character- 
ized by a Coulomb repulsion that is comparable to or 
larger than hopping constants. One of the simplest mod- 
els for the description of such systems is the Hubbard 
model-'^'^ with the Hamiltonian 



H 



u 



(1) 



. In the case of 
nmi, it is reasonable 



where tnm is the hopping constants, the operator ajj^. 
creates an electron on the site n with the spin projec- 
tion tr = ±1, J7 is the on-site Coulomb repulsion, the 
electron number operator Una = ajj^, 
strong electron correlations, U ^ \t 
to use a perturbation expansion in powers of the hop- 
ping constants for investigating Hamiltonian ^ . Appar- 
ently the first such expansion was considered in Ref. jjj. 
The fu rthe r development of this approach was given in 
Refs. 000110 where the diagram technique for Hubbard 
operators was developed and used for investigating the 
Mott transition, the magnetic phase diagram and the su- 
perconducting transition in the Hubbard model. 

The rules of the diagram technique for Hubbard oper- 
ators are rather intricate. Besides, these rules and the 
graphical representation of the expansion vary depend- 
ing on the choice of the operator precedence. The dia- 
gram technique suggested in Refs. 10 11 12 13 is free from 
these defects. In this approach the power expansion for 
Green's function of electron operators ana and aj,^ is con- 
sidered and the terms of the expansion are expressed as 
site cumulants of these operators. Based on this diagram 
technique the equations of the Larkin type^** for Green's 
function were derivediiSiiSii^ 

However, the application of this approach runs into 
problems. In particular, at half-filling the spectral weight 
appears to be negative near frequencies ivd = ± This 
drawback is connected with divergencies in cumulants at 
these frequencies. As can be seen from formulas given 
below, all higher-order cumulants have such divergencies 
at LOd with sign-changing residues, which are expected 



to compensate each other in the entire series. If, as in 
Ref. UJ, only some subset of diagrams is taken into ac- 
count the divergencies of different orders can be com- 
pensated only accidentally. On the other hand, at fre- 
quencies neighboring to cud cumulants are regular. If a 
selected subset of diagrams is expected to give a cor- 
rect estimate of the entire series for these frequencies the 
values at ujd can be corrected by an interpolation using 
results for the regular regions. This procedure is applied 
in the present work where the one-loop approximation is 
used ~ as irreducible diagrams in Larkin's equation those 
appearing in the first two orders of the perturbation ex- 
pansion are adopted. One of these diagrams contains a 
loop formed by a hopping line. It is shown that for the 
case of a square lattice, nearest-neighbor hopping and 
half-filling the electron spectrum consists of four bands. 
These band structure and the calculated shapes of the 
electron spectral function are close to those obtained in 
Monte-Carlo calculation9iS*iS*ii and recently in the clus- 
ter perturbationi^ and the two-particle self-consistenU^ 
calculations. It is shown also that the maxima forming 
the bands are of a dissimilar origin in different regions of 
the Brillouin zone. 

In the following section the perturbation expansion for 
the electron Green's function is formulated in the form 
convenient for calculations and Larkin's equation is de- 
rived. In Sec. Ill the equations of the previous section 
are used for calculating the spectral function and the 
obtained results are discussed. Concluding remarks are 
presented in Sec. IV. 



II. DIAGRAM TECHNIQUE 

We consider Green's function 

G(nr,n'T')- (Ta„,,(r')an,(r)), 



(2) 



where the angular brackets denote the statistical averag- 



ing with the Hamiltonian Ti — H — ^ 



IX is the 



chemical potential, T is the time-ordering operator which 
arranges other operators from right to left in ascending 



2 



order of times r, ano-(T) = exp('WT)a„cr exp(— Hr) and 
anaiT) — exp{TiT)al^^exp{—TiT). Choosing 



where 



U 



(3) 



as the unperturbed Hamiltonian and the perturbation, 
respectively, and using the known^-5. expansion for the 
evolution operator we get 



G(n'r',nr) = ^ 

A;=0 



A:! 



dri . . . dTk 



E 



nin'^CTi nfcn'j.(Tfc 

X (Ta„'o-(T')a„o.(T)a„'^cri(Ti)anio-i(Ti) • ■ • 
X a„;^<T,(T-fc)a„,^,(Tfc))oc, (4) 

where f3 — 1/T is the inverse temperature, the subscript 
"0" near the angular bracket indicates that the averaging 
and time dependencies of operators are determined with 
the Hamiltonian Hq. The subscript "c" indicates that 
terms which split into two or more disconnected averages 
have to be dropped out. 

The Hamiltonian Hq is diagonal in the site represen- 
tation. Therefore the average in the right-hand side of 
Eq. Q splits into averages belonging to separate lattice 
sites. To be nonzero these latter averages have to contain 
equal number of creation and annihilation operators. Let 
us consider the average which appears in the first order: 
(Tan'o-(T')a„o-(T)a„'^cri(Ti)a„io-i(Ti)) (for short here and 
below subscripts "0" and "c" are omitted). For this av- 
erage to be nonvanishing operators have to belong either 
to the same site or to two different sites, 

(Ta„'o-(T')a„o-(T)a„'^cri (r{ )a„jo-i (n)) = 

{Tdna{T')ana{T) {t[ ) anai (n )) f^nn' 6nn[ Snn^ 

+ {Tancr{T')ana{T)) {Tdmai {T[)an^^^ (ti)) 

X(5nn"5ninj (1 '^nni ) 
-(Tanio-(T')anio-i (ti)) (TUnai (T{)a„cr(T)) 

The multiplier 1— ^nni ensures that n rii in the two last 
terms (the case n = ni is taken into account by the first 
term in the right-hand side) . After the rearrangement of 
the terms we find 

{Tan'cr{T')ancriT)an'^ai (TOaniCTi (n)) = 
K2{t' 

+Ki(T'a, rcr)Ki(T{a-i,Ti(Ti)5„„/(5„j„'^ 

-Ki {t'cJ, Ti(Ti)Ki (t{ CTi , Ta)5n'm ^nni , 



K^{T'o',Ta) = {Ta„{T')a,{T))5„ 
K2{T'a,Ta,T[ai,Tiai) = 

{Tda{T')a^{T)dai (rOflcri (ti)) 
-Kiir'a, Ta)Ki{T[ai,Tiai) 
+Ki{t' (j,Ti(7i)Ki(T[ai,Ta) 



(5) 



are cumulantsSi of the first and second order, respec- 
tively. All operators in cumulants belong to the same 
lattice site. Due to the translational symmetry of Hq in 
Eq. ||2Jl the cumulants do not depend on the site index 
which is therefore omitted in Eq. 

Averages which appear in higher-order terms of Eq. Q 
can be transformed in the same manner. The average in 
the fc-th order term which contains fc -I- 1 creation and 
annihilation operators is represented by the sum of all 
possible products of cumulants with the sum of orders 
equal to fc-l- 1. All possible distributions of operators be- 
tween the cumulants in these products have to be taken 
into account. The sign of a term in this sum is deter- 
mined by the number of permutations of fermion opera- 
tors, which bring the sequence of operators in the initial 
average to that in the term. Actually the above state- 
ments determine rules of the diagram technique. Ad- 
ditionally we have to take into account the presence of 
topologically equivalent terms ~ terms which differ only 
by permutation of operators Hiiji) in Eq. Q. Since 
these terms are equal, in the expansion only one of them 
can be taken into account with the prefactor v = j /kl 
where j is the number of topologically equivalent terms. 
Following Ref. |ll| in diagrams we denote a cumulant by 
a circle and the hopping constant inn' by a line directed 
from n' to n. The external operators dn'^iT') and ancr(r) 
are denoted by directed lines leaving from and entering 
into a cumulant. The order of a cumulant is equal to a 
number of incoming or outgoing lines. Summations and 
integrations over the internal indices n^, n[, cr,; and 
are carried out. Since site indices of operators included 
in a cumulant coincide, some site summations disappear. 
Also some summations over ai get lost, because in any 
cumulant spin indices of creation and annihilation oper- 
ators have to match. Taking into account the multiplier 
(—1)*^ in Eq. the sign of the diagram is equal to 
(— f )' where I is the number of loops formed by hopping 
lines. Figure n demonstrates connected diagrams of the 
first four orders of the power expansion ^ with their 
signs and prefactors. Here the thick line with arrow in 
the left-hand side of the equation is the total Green's 
function. Notice that if we set t„n — 0, contributions of 
the diagrams (b), (d)-(f), (i)-(n), (p), (t), and (u) van- 
ish. However, below a rcnormalizcd hopping parameter 
will be introduced which is nonzero for coinciding site in- 
dices and therefore the mentioned diagrams are retained 
in Fig.n 

All diagrams can be separated into two categories - 
those which can be divided into two parts by cutting 
some hopping line and those which cannot be divided in 
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transformation 



(a) (b) (c) 



(d) 



(f) (g) (h) 



G(k, it^i) = V e*''^"'-"' / e*"'(^'-^)G(n'T', nr), 

where = {21 + l)7rT is the Matsubara frequency, sim- 
pUfies significantly these equations: 




K{k, tiui) 



1 - <kA'(k, iuji) 



(8) 



(9) 



(q) (r) 



(u) (V) 

FIG. 1: Diagrams of the first four orders of expansion (|1J. 

this way.^-^'' These latter diagrams are referred to as ir- 
reducible diagrams. In Fig. ^the diagrams (c), (e), (f), 
(h), (j), (k), (n), (p), and (r)-(v) belong to the first cate- 
gory, while the others are from the second category. The 
former diagrams can be constructed from the irreducible 
diagrams by connecting them with hopping lines. No- 
tice that a prefactor v of some composite diagram is the 
product of prefactors of its irreducible parts [cf., e.g., the 
irreducible diagrams (a), (d) and the composed diagrams 
(j), (k) in Fig.^. If we denote the sum of all irreducible 
diagrams as K{n'T' ^nr), the equation for Green's func- 
tion can be written as 



In the approximation used below the total collection of 
irreducible diagrams K{\i,iu!i) is substituted by the sum 
of the two simplest diagrams (a) and (b) appearing in the 
first two orders of the perturbation theory. Due to the 
form of the latter diagram this approximation is referred 
to as the one-loop approximation. In the diagram (b) the 
hopping line is renormalized in accordance with Eq. 
Thus, 



K{iuji) = Ki{iuJi) -T'^K2{iuJia,iuJi^ai,iuJi^ai) 



1 



k 



(10) 



where Ki{iu!i) and 



dr' drdT'idTie 



X K2 (t'ct, TCr, t[ 0-1 , Tl (Tl ) 



G(n'T',nT) = K{n'T',nT)+ ^ / dTiK{n'T' ,niTi) are the Fourier transforms of cumulants jSJl, N is the 

mn; number of sites and we set X^k^k — 0. Notice that in 



(6) 



This equation has the form of the Larkin equation»i4 

The partial summation can be carried out in the hop- 
ping lines of cumulants by inserting the irreducible dia- 
grams into these lines. In doing so the hopping constant 
tnn' in the respective formulas is substituted by 

e(nT,nV) = t„„'(5(T-T') 



f 

V inn; / dTiK{n[T,niTi) 
, Jo 



number of sites and we set X^k^k — 0. Notice that 
this approximation K does not depend on momentum. 

Now we need to calculate the cumulants in Eq. H1U|I . 
To do this it is convenient to introduce the Hubbard op- 
erators X^^ = |in)(jn| where |m) are eigenvectors of site 
Hamiltonians forming Hq, Eq. ||2Jl. For each site there 
are four such states: the empty state |0n) with the en- 
ergy Eq = 0, the two degenerate singly occupied states 
\an) with the energy Ei ~ — /i and the doubly occupied 
state |2n) with the energy E2 = U — 2/i. The Hubbard 
operators are connected by the relations 



iL = XZ° + aXl'-" 



(11) 



X 9(niri, n'r'). 



(7) 



For the diagram (b) in Fig. ^ inserting irreducible dia- 
grams into the hopping line leads to the diagrams (g), 
(m), and (q). 

Due to the translation and time invariance of the prob- 
lem the quantities in Eqs. 10 and Q depend only on dif- 
ferences of site indices and times. The use of the Fourier 



with the creation and annihilation operators. The com- 
mutation relations for the Hubbard operators are easily 
derived from their definition. Using Eq. (|ll|l the first 
cumulant in Eq. can be computed straightforwardly: 



Ki{iuji) ^ — 
^0 



1 (e-l"^' 



iuJi - EcrO 



iuJi - E2a 



(12) 
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2g-^s„ _|_ g-)3E2 ^jjg partition where gij^icui) = (iuji + Eij)~^ is the Fourier transform 

" of functions (HH) and lfT3|) . 

In this work the case of half-fiUing is considered. In this 
case /i = -J and cumulants (|12|1 and are significantly 
simplified if we additionally suppose that T <C C/ : 



where Zq = er/^^" 

function, Eij — Ei — Ej. As indicated in Refs. ISl llOlll,: 
if K^k, iuJi) is approximated by this cumulant the result 
corresponds to the Hubbard-I approximation.^ 

To calculate K2 it is convenient to use Wick's theorem 
for Hubbard operators (^iSiS 

{TXa,, (n) . . . X,, (t,)X,(t)X„,^, (t,+i) . . . X„„ (r„)) 
■(-l)^^g„(r-Tfe)(rX„,(ri)... 



fc=i 



(13) 



where the averaging and time dependencies of operators 
are determined by the Hamiltonian Hq, a is the index 
combining the state and site indices of the Hubbard op- 
erator. If Xa is a fermion operator (X^*^, X'^^ and their 
conjugates), Pk is the number of permutation with other 
fermion operators which is necessary to transfer the oper- 
ator Xa from its position in the left-hand side of Eq. H13|) 
to the position in the right-hand side. In this case 

PT|e/5^S r<0, 

where i and j are the state indices of Xa- If Xa is a 
boson operator (X™, X^^, X'^'^', X"^, and X^O), Pj, = 
and 

gg-.- f 1, r >0, 



1 



In Eq. 113|l . [Xa^ , ^a]± denotes an anticommutator when 
both operators are of fermion type and a commutator in 
other cases. 

The substitution of Eq. ifTT)) in (T dcr{T')aa{T)aai{T{) 
xaCTj(Ti)) in K2, Eq. |(SJ, leads to six nonvanishing aver- 
ages of Hubbard operators (such averages are nonzero if 
the numbers of the operators X'^°' and X°''^ coincide in 
them, and the same is true for the pair and X'^°'). 
Applying Wick's theorem (|13() to these averages the num- 
ber of operators in them is sequentially decreased until 
only time-independent operators are retained. For Hq in 
Eq. © these are X°°, X""' , and X^^. Their averages 
are easily calculated. As a result, after some algebra we 
find 

'^K2{iL0ia,iLOi^CFi,U}i^ai) = -Z^'^U {e~'^^'> g(i^{iLOi) 



-I3E, 



y.[ga2{it^i) + gc,2{iuji^)\ +e '^^^[gQc,{it^i)ga2{it^i) 
x{gaa{ii^h) - ga2{iuJh))'^ + gtia{i^h)ga2{iuJi^) 
X {gU^u:,) + gl^oJi))] ] Z,'U^(35u, ^,-Ke.+e.) 
^2e-/5(Bo+Bi) + 3e-2/5£i ^ 2e-^'^'''+''-^)gl„{iL0i) 



\gQa{it^l)ga2{it^l)gQa{it^h)ga2iil^h), 



(16) 



Ki{iLUi) 



lUJl 



U 



2 V[(-o^-(#) ] 
[(z.o^-(#)Y[(*c..j2-(l)^ 



(17) 



-f3Sih 



The first term in the equation for K2 is the even 
function of iwi^. It does not contribute to the sum over 
li in Eq. (fTU)) . since X^k ^k'-^C*^' ) supposed to be an 
odd function of ioJi^. Thus, finally Eq. acquires the 
form 



lUJl 



+ 



^(ill-^lE4G(k,..,).(18) 



Analogous equations were obtained in Refs. I8liril3l (in 
Ref. 1^ the last term in the right-hand side of the equa- 
tion is three times smaller, because diagrams with spin 
propagators were not taken into account there). 

Let us turn to real frequencies by substituting ioji with 
z = uj -\- irj where 77 is a small positive constant which 
affords an artificial broadening. Results given in the next 
section were calculated with G(k, z) in the right-hand 
side of Eq. \i'6\ taken from the Hubbard-I approximation. 
As mentioned, this Green's function is obtained if K^kuj) 
in Eq. © is approximated by the first cumulant (|12ll . For 
half-filling we find 



G(k,z)^l 1+ ^ 

2 I .JWTti ^-ei,k 



1 




ei.k = ^ ( tk + \/C/2 + tl,, 



(19) 



Below the two-dimensional square lattice is considered. 
It is supposed that only the hopping constants between 
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(0,Ji)k 



FIG. 2: The imaginary part of K{u}) calculated using 



Eqs. (EHJ and ^ for a 100x100 lattice, t 



i U and 



rj = 0.01(7 (the dashed line). The solid line shows the cor- 
rected ^K{ijo) (see text). 



nearest neighbor sites t are nonzero which gives ik — 
2t[cos{kx) -|-cos(fcj,)] where the intersite distance is taken 
as the unit of length. 




III. SPECTRAL FUNCTION 

Figure 12 demonstrates '^K{lo) calculated with the use 
of Eqs. H18|) and H19I) . The change to real frequencies car- 
ried out in the previous section converts the Matsubara 
function Q into the retarded Green's functioni^ It is an 
analytic function in the upper half-plane which requires 
that '^K{uo) be negative. As seen from Fig.|21 this condi- 
tion is violated at w = This difficulty of the consid- 
ered approximation was indicated in Ref . |3 The prob- 
lem is connected with divergencies at = A* ^-nd ^ — U 
introduced by functions goa-^Lu) and ga-2{i-u) in the above 
formulas. As can be seen from the procedure of calculat- 
ing the cumulants in the previous section, these functions 
and divergencies with sign-changing residues will appear 
in all orders of the perturbation expansion Q . It can be 
expected that in the entire series the divergencies com- 
pensate each other, however, in the considered subset of 
terms such compensation does not occur. Nevertheless, 
as seen from Fig. |21 at frequencies neighboring to uid cu- 
mulants are regular and, if the used subset of diagrams 
is expected to give a correct estimate of the entire series 
for these frequencies, the value of ^K{uj) at the singular 
frequencies can be reconstructed using an interpolation 
and its values in the regular region. An example of such 
interpolation is given in Fig. |21 

The function K{z) has to be analytic in the upper half- 
plane also and therefore its real part can be calculated 
from its imaginary part using the Kramers-Kronig rela- 
tions. We use this way with the interpolated 3i^(w) to 
avoid the influence of the divergencies on '3iK{u!). How- 



FIG. 3; The spectral function A(kLu) calculated for momenta 
along the route shown in the inset. An 8x8 lattice, t — ~jU 
and ri = 0.2)7. 

ever, the use of the interpolation overrates somewhat val- 
ues of |3i4r(a;)| which leads to the overestimation of the 
tails in the real part. To correct this defect the interpo- 
lated K{u!) is scaled so that its real part for frequencies 
Iti-"! ^ coincide with the values obtained from Eq. I|18|l . 
It is worth noting that such obtained spectral function 

A{kLj) = ^-^G{kuj) 

= A (20) 

satisfies the normalization condition dujA{]iuj) = 1 
with an accuracy better than 0.01. 

The spectral function obtained in this way with the 
use of Eqs. ((TH|l . ifT^ . and for momenta in the Bril- 
louin zone of an 8x8 lattice is shown in Fig.O In these 
calculations the small lattice and a comparatively large 
artificial broadening were used to reproduce the condi- 
tions of Monte-Carlo calculations carried out in Ref. 12^ . 
Comparing Fig. with Fig. 5 in that work we find close 
similarity of these two groups of spectra - not only the 
general shapes of spectra are close but also the position 
of maxima in them. Some discrepancies can be also no- 
ticed. In particular, in the spectrum for k = (^, j) in 
Fig. 131 the intensity of maximum at a positive fre que ncy 
is larger than in the respective spectrum in Ref. |23 (in 
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this latter spectrum this maximum looks like a barely 
perceptible peculiarity). 

The spectral function for a larger U is shown in Fig.0] 
where the momenta along the symmetry lines of the 
square Brillouin zone were chosen. Let us compare these 
results with the spectra in Fig. 1 in Ref. |0 which were 
obtained in Monte-Carlo calculations with the same ra- 
tio U/\t\ in an 8x8 lattice. We found that the results 
shown in Fig. 01 change only slightly with decreasing the 
lattice size to this value. The comparison has to be car- 
ried out with low-temperature spectra of Ref. llTl since 
Eq. H17|l used in our calculations was derived for the case 
T <C C/. Notice that in this limit at half-filling the spec- 
tral function in the used approximation does not depend 
on temperature. Again we find that our calculated spec- 
tra are close to those obtained in the Monte-Carlo cal- 
culations, though there are some differences in relative 
intensity of maxima. In our spectra some maxima re- 
solved in the Monte-Carlo calculations look like shoulders 
of main maxima and a decrease of the artificial broaden- 
ing T] does not convert the shoulders to maxima. On 
the other hand, the maxima at w w 2\t\ in the spectra 
for momenta near (0,0) and aX u ~ — 2|i| in the spectra 
for momenta near (tt, tt) which are well resolved in our 
results are barely perceptible in Fig. 1 of Ref. ^3- The 
possible reason for these differences can be established 
from the comparison of Figs. |31 and ^ The respective 
maxima in the former figure for the smaller U are less in- 
tensive. Apparently, like the Hubbard-I approximation, 
the used approach somewhat overestimates the electron 
correlations. Nevertheless, Fig. ^ reproduces correctly 
the most essential features of the spectra obtained by 
the Monte-Carlo method. 

One of these features recently reproduced also in 
the cluster perturbationiSi and the two-particle self- 
consistent^^ calculations is the four-band structure of the 
spectrum. In Fig.^these bands are located near frequen- 
cies Lo = ±3|t| and oj — ±6|t|. The high-frequency bands 
are observed near the (0,0) and (tTjTt) points, while the 
low-frequency bands have a larger existence domain in 
our calculations. Similar four bands can be identifies in 
Fig. 131 A more detailed treatment of the low-frequency 
bands shows that either of them in turn contain two com- 
ponents. Maxima in A{kijj) arise by two reasons - due 
to the smallness of the denominator in Eq. H20|) when a 
frequency satisfying the equation 1 — t\JkK{ijj) = falls 
into a region of a small damping |3_ftr(a;)| and due to 
the strong frequency dependence of 5-ftr(w) in the nu- 
merator of Eq. (|20|l . The high-frequency bands and the 
low-frequency bands in the central part and at the pe- 
riphery of the Brillouin zone arise due to the first reason. 
As seen from Fig. 12 the regions of small damping are lo- 
cated on either side of the two maxima in '^Kiuj) which 
gives rise to the well-separated four bands. However, for 
momenta which are close to the boundary of the magnetic 
Brillouin zone the equation 1 — t\JkK{ijj) — has no so- 
lutions due to a small or vanishing tk- In this case the 
two maxima of Ai^kw) originate from the two maxima in 




-12 -8 -4 4 8 12 



ca/\t\ 



T r 




(0,0) (ji,o) (7i,7i) (0,0) 



FIG. 4: The upper panel: the spectral function A(ka;) cal- 
culated for momenta along the symmetry lines of the square 
Brillouin zone. The lower panel: the dispersion of maxima 
in the upper panel shown as a gray-scale plot with darker ar- 
eas corresponding to larger intensity of maxima. A 100x100 
lattice, t = —jU and rj = 0.117. 
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1 3/^(0-1)1. They are located in the same frequency ranges 
as the low-frequency maxima originated from the reso- 
nance denominator in Eq. (|2()(l . Therefore they can be 
combined in a common dispersion band, although widths 
of the maxima of different origin can significantly vary for 
a smaller artificial broadening. 

As can be seen from Eqs. (jS)), (|18|) . and H19|) . the used 
approximation does not describe the Mott transition - 
for 7/ — > the density of states = A^"^ Sk^l^'^) 
at w = vanishes for any finite U . This flaw can be 
remedied by going to the next approximation. The use 
of approximation H19|l in Eq. (|18|) means that only dia- 
grams (a) in Fig. ^ arc inserted in the hopping line of the 
diagram (b). If these latter diagrams are also inserted 
in the line, Eq. (|18|l is transformed to the following self- 
consistent equation for the irreducible part: 



+ [,2_(a)YiV^l-<kifW ^^^^ 

For the initial electron dispersion described by the 
semielliptical density of states 

I 0, |c.|>i3, 

where B is the halfwidth of the band, Eq. (|21(l can be 
solved analytically and the condition for the disappear- 
ance of the gap at w = in piiS) can be found. It happens 
when U < Uc — V^B and the gap opens when U > Uc- 
Analogous results for somewhat different equat ions and 
initial dispersions were obtained in Refs. rTl lSllSi 



parameters. The expansion in powers of these latter pa- 
rameters can be expressed in terms of cumulants in the 
same manner as discussed above. Now there are distinct 
cumulants which belong to different site states charac- 
terized by dissimilar parameters of the repulsion and the 
level energy. These cumulants are described by formu- 
las similar to Eqs. H12|l and H16|) . For example, for the 
Emery model^^ which describes oxygen 2pa- and copper 
3dj;2_y2 orbitals of Cu-0 planes in high-Tc superconduc- 
tors there are two types of cumulants corresponding to 
these states. Diagrams of the lowest orders, e.g., for 
Green's function on copper sites resemble those shown 
in Fig. ^ where "oxygen" cumulants are included in hop- 
ping lines. Equations of the type of Eq. and partial 
summations similar to Eq. (PJ can be derived also in this 
case. 

In summary, the diagram technique for the one-band 
Hubbard model was formulated for the case of moderate 
to strong Hubbard repulsion. The expansion in pow- 
ers of the hopping constant is expressed in terms of site 
cumulants of electron creation and annihilation opera- 
tors. For Green's function the equation of the Larkin 
type was derived and solved for the case of two dimen- 
sions, nearest-neighbor hopping and half-filling. The ob- 
tained spectrum contains four bands. It was shown that 
maxima of the spectral function forming the bands are 
of dissimilar origin in different regions of the Brillouin 
zone - in its central part and at the periphery the max- 
ima stem from the resonance denominator in the formula 
for the spectral function, while near the boundary of the 
magnetic Brillouin zone they originate from the strong 
frequency dependence of the irreducible part in the nu- 
merator of this formula. The obtained band structure 
and the shapes of the spectral function are close to those 
found in the Monte Carlo calculations. 



IV. CONCLUSION 



The considered diagram technique is very promising for 
a generalization to many-band Hubbard models for which 
energy parameters of the one-site parts of the Hamilto- 
nians exceed or at least are comparable to the intersite 
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